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ABSTRACT 

We describe a major upgrade of a Monte Carlo code which has previously been used for 
many studies of dense star clusters. We outline the steps needed in order to calibrate the 
results of the new Monte Carlo code against iV-body simulations for large systems, up to 
TV = 200000. The new version of the Monte Carlo code (called MOCCA), in addition to the 
old version, incorporates direct FewBody inte grator for three- and four-bod y interactions, and 
new treatment of the escape process based on lFokushige & Heggid (|2000|) . Now stars which 
fulfil the escape criterion are not removed immediately, but can stay in the system for a certain 
time which depends on the excess of the energy of a star above the critical energy. They are 
called potential escapers. FewBody integrator allows to follow all interaction channels, which 
are important for the rate of creation of various types of objects observed in star clusters, and 
assures that the energy generation by binaries is treated in a meaner similar to the N-body 
model. 

There are at most three parameters which have to be adjusted against N-body simula- 
tions for large A^. Two (or one, depends on the chosen approach) connected with the escape 
process and one responsible for determination of the interaction probabilities. The adopted 
free parameters are independent on A^. They allow MOCCA code to reproduce N-body re- 
sults, in a reasonably precision, not only for the rate of cluster evolution and the cluster mass 
distribution, but also for the detailed distributions of mass and binding energy of binaries. 
Additionally, the code can follow the rate of formation of blue stragglers and black hole - 
black hole binaries. The code computes interactions between binaries and single stars up to a 
maximum separation rpmax, and it is found that the MOCCA code needs rather large value of 
Tpmax to get agreement with the N-body simulations. 

The MOCCA code is at present the most advanced code for simulations of real star clus- 
ters. It can follow the cluster evolution in details comparable to N-body code, but orders of 
magnitude faster 
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1 INTRODUCTION 

This is the second paper in a new series of papers in which we at- 
tempt to describe the development of the MOCCA (MOnte Carlo 
Cluster simulAtor) code and its applicatio n to the simulations of 
star cluster evolution. The first in the series jHvDki & Gierszl201lh 
described in detail rece nt developments of the previous version of 
the Monte Carlo code ( iGiersz. Heggie. & HurlevI ll2008l). and ref- 
erences therein) and the first results of simulations concerning blue 
stragglers (BSS) in an evolving star cluster environment. In this pa- 
per, we further develop the code and perform a very detailed com- 



* E-mail: mig@camk.edu.pl (MG) 
© 2002 RAS 



parison with N-body simulations of large A'' stellar systems up to 

TV = 2 X 10^ 



The MOCCA code l lHvpki & Gie"rs3l2oT J) is at present the 
most advanced numerical code for stellar dynamical simulations, 
capable to follow evolution of real star clusters in detail com- 
parable to N-body simulations, but orders of magnitude faster 
(several hours for N = 2 x 10®). The dynamical ingredients 
of the Mo nte Car l o code are essentially the same as those de - 
scribed in lOiersj ( l2006h and lOiersz. Heggie. & HurlevI (l2008h, 
whose code embodies several features introduced by Stodolkiewiczl 
( Il986l). wh ose co de was in turn based on that originally devised 
bv iHenonI ( Il97ll) . Two main features distinguish the MOCCA 
code from the previous version of the Monte Carlo code: (i) it 
now incorporates dynamical interactions between binary and sin- 
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gle stars and between pairs of binari e s bas ed on the FewBody 
integrator developed by IFregeau et all ( l2004i) : (ii) it replaces the 
treatment of the e scape process in the static tidal field based on 



Baumgardt! (120011) b y one i n accordance the theory proposed by 



Fokushige & Heggid Jioooh . The escape process is not instanta- 



neous any more, an object needs time to find its way around the 
Lagrangian points Li and L2 to escape. The MOCCA code in- 
corporates most of the processes which are important during stel- 
lar system evolution e.g.: the relaxation process, the main en- 
gi ne of dynamical clus t er evo lution; stellar evolution according 
to iHurley. Pols. & Tou 3 1 I2OOOI) f or the evolution of s i ngle s tars , 
supplemented by the methods of IHurlev. Tout. & Pols! j2002h for 
internal evolution of binary stars (BSE code) and also simple 
approac h for colliding stars accordi ng to McScatter interface to 
BSE bv lHeggie. Portegies Zwart. & H urley (200^); the escape pro- 
cess in the static tidal field of the parent galaxy; direct few 
body integrator to follow interactions between binaries and single 
stars and other bin aries; mass segregated initial cluster con figu- 
ration acc ording to Baumgardt, De Marchi. & Kroupj 1 I2OO8I) and 
ISubr. Kroup a & Baumgardt ( 2008). 

There are several factors which motivate this work. Star clus- 
ters are the focus of many intensive observational campaigns 
(e.g. iBedin et alj I2001I; iBedin. Piotto et al. 1 12003|; IGrindlay et al.l 



200 li; 



Piotto et al 



Richer et al 



20041 ; 



2002; 'Kaliari et al. 2003; Kafka et a^ 
Anderson et al., 2006 ; Milone et al.. 20 11 



liay et ai.l 
1IJI 2OO4 



and 

references therein), which are now turning to the examination of the 
parameters of their populations of different kinds of binaries, BSS 
and other "peculiar" objects. Dynamical models are needed for the 
design and interpretation of observational programmes: how is the 
period distribution and the spatial distribution of binaries affected 
by dynamical evolution? What is the influence of environment and 
dynamical evolution on the formation of "peculiar" objects? Un- 
derstanding the abundance, spatial distribution and channels of for- 
mation of BSS can only be attempted by a technique which follows 
simultaneously both BSS dynamics and internal evolution. While 
A'^-body technique may ultimately be the method of choice for such 
studies, systems of the size of most globular clusters are likely to 
remain beyond reach for some years, simply because of the num- 
ber of stars and the size of the binary populatio n. After all, it is only 
recently that the "hardest" open clusters M67 jHurley et alj|200a) 
and "easiest", loosely b ound and distant globular cluster Palomar 
14 jZonoozi et alj201ll) have been modelled at the necessary level 
of sophistication. To efficiently compute detailed models of large 
star clusters and to investigate the influence of initial parameters on 
a cluster's global and local observational properties we need a tech- 
nique which is much faster than the N-body code and at the same 
time can give the same level of information about every object in 
the cluster as the N-body code does. The MOCCA code is such a 
technique. 

One of the drawbacks of the non direct techniques (also the 
Monte Carlo one) compared to the N-body model is a necessity 
of use of free parameters which try to describe the complexity of 
physical processes naturally covered in the direct code. The most 
important free parameters (from the point of view of MOCCA) are 
connected with the relaxation process (7 coefficient in the Coulomb 
logarithm), escape process in the static tidal field and dynamical in- 
teractions between different objects (where the parameter, rpmax, 
is the maximum pericentre distance between interacting objects for 
which few-body interactions are calculated explicitly). The usual 
method to determine the free parameters is a comparison with 
the results of N-body simulations. For the previous version of the 
Monte Carlo code the comparison was done only for small TV sys- 



tems (up to = 24000). The code was successfully used to simu- 
late e voluti on of real star clusters: M 67 dOiersz. Heggie. & HurlevI 
2008 1) M4 llHeggie & Gier"s^l2008al) . NGC6397 taiersz & Heggie 



2009l ; lHeggie & Gierszll200a) and 47Tuc dOiersz & Heggieil201ll) 



Despite those successes there were some doubts connected with 
the A'' s caling of the escape process, which implementation was 
based on lBaumgard] faoii) . To fully trust the MOCCA code it has 
to be tested for larger N, and not only for the global parameters 
like evolution of the total cluster mass or Lagrangian radii, but also 
against the properties and spatial distributions of binaries and BSS. 
That kind of comparison will show how far we can trust results of 
Monte Carlo simulations and which processes cannot be properly 
described in the framework of the MOCCA code. 

This paper begins in Sec. 2 with a summary of the features 
which have been added to the Monte Carlo scheme during the con- 
struction of the new version of the code, the MOCCA code. We also 
show there how we calibrate the free parameters of the MOCCA 
code with results of A^-body simulations. Next (Sec. 3) we describe 
the similarities and differences between MOCCA and N-body sim- 
ulation and discuss the possible reasons for that. The final section 
summarises our conclusions, and discusses some limitations and 
future developments of the MOCCA code. 



2 TECHNIQUE 

The MOCCA code l lHypki & Gierszl I2OIII) is an updated ver- 
sion of the Monte Carlo code develope d in iGierszi ( 1 19981 . 1200 iL 
I2OO6I) ; iGiersz. Heggie. & HurlevI 1 I2OO8I) . In addition to the de- 
scription of the relaxation process which is responsible for the 
dynamical system evolution it includes synthetic stellar evolu- 
tion of single and b i nary s tars using prescriptions d e scribe d by 
IHurlev. Pols. & TouJ ( I2OO0I) and IHurlev. Tout. & Polsl ( I2OO2I) and 
direct integration pr ocedures f or small A*' subsystems based on the 
FewBody code iFregeau et ^ 12004) . One of the more important 
updates is a better descript ion of the escape process according to 
jFokushige & HeggidbOOd) . Now the escape of an object from the 
sy stem is not instantaneous any more, but time delayed. The theory 
of iFokushige & Heggid ( I2OO0I) incorporates a number of parame- 
ters which they had to determine empirically, and which depend 
on the system under consideration; here we shall determine these 
parameters by comparing the results with those of N-body simu- 
lations. In the N-body code formation of binaries and subsequent 
their interactions follow naturally from the movement of stars in 
the system under the mutual gravitational force. In the MOCCA 
code first we need to check if the interaction is due be computing 
its probability and then if it is so we execute the direct integration 
procedure for small A*" (3- or 4-body) subsystem to find out the 
outcome of the interaction. The interaction probability depends, 
among others, on the maximum value of the pericentre distance, 
Tpmax- The larger this distance the larger the number of interac- 
tions, which are weaker on average. Choosing a proper value of 
rpmax is crucial for a balance between code efficiency and its ac- 
curacy, e.g. the number of BSS observed in the system strongly 
depends on rpmax- 



2.1 Delayed escape 

As it was poin ted out in IFokushige & Heggid 1 I2OOOI) and 
IBaumgardtl ( 1200 ll) the process of escape from a cluster in a steady 
tidal field is extremely complicated. Some stars which fulfil the en- 
ergy criterion for escape (binding energy of the star greater than 
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Table 1. Initial conditions for N-body simulations 



Cluster 


N=24000 (M67) 


N= 100000 (NGC6397) 


N=200000 


Ns 


12000 


95000 


195000 


Nb 


12000 


5000 


5000 


Binary fraction 


0.5 


0.05 


0.025 


M(0) 


1.869 X IO-'Mq 


5.177 X lO'' Mq 


1.001 X lO^Mt 


Initial model 


Plummer 


Plummer 


Plummer 


Initial tidal radius 


32.2pc 


52.4pc 


35.8pc 


IMF of stars 


Kroupa" 


Kroupa" 


Kroupa*^ 


IMF of binaiies 


Kroupa** 


Kroupa*" 


Kroupa'' 


Mass ratio 


Uniform 


Uniform 


Uniform 


Binary eccentricities 


thermal'^ 


thermal 


thermal 


Binary semi-major axes 


Uniform'' 


Uniform'' 


Uniform'' 


SN kick distribution 


Gaussian*^ 


Gaussian"^ 


flat J' 


Metallicity 


0.02 


0.001 


0.001 



" iKroupa. Tout & Gihnon 
'' iKioupa. Grhnore & Toui 



199: 



with mass range between 0.1 and 5OM0 

19911) eq.(l) wit h mass range between 0.2 and IOOMq 

'^ thermal distribution modified according to lHurlev et all )2005h eq.(l) 

'' uniformly distributed in the logarithm in the range 2(iJi -|- R2) to 50AU 

Gaussian distribution with a = 190 km/s 
f uniform distribution with kick velocities between to 100 km/s 



the critical energy Ecru = —1.5{GM /rt), where G is the gravi- 
tational consta nt, M is the total mass and rt is the tidal radius, see 
ISDitzeil l fl987l) ) can still be trapped inside the potential well. Some 
of those stars can be scattered back to lower energy before they es- 
cape from the system. These two factors cause the cluster lifetime 
to scale nonlinearly with relaxation time for tidally limited clus- 
ters (Baumgardt 2001;), in contrast with what would be expected 
from the standard theory. The efficiency of this effect decreases as 
the number of stars increases. To account for the process described 
above in the previous version of the Monte Carlo code an addi- 
tional free parameter a was introduced ( iGiersz. Heggie. & Hurley! 
I2OO8I) . The critical energy for escaping stars was approximated by: 
E,„t = -xud{GM/rt), where xua = 1.5 - a{ln{yN)/N)^^\ 
a was approximated by 2.5 and 7 is the coefficient in the Coulomb 

logarithm equal to 0.11 for equal mass case and 0.02 for un- 

I : ' ' 1 

equal mass case (see IGiersz. Heggie. & Hurley 2008. a nd refer- 
ence t herein). This prescription was also tested by Chat teriee et al.l 
( I2OIOI) in their Monte Carlo simulations of star clusters. They are 
three drawbacks of this approach: (i) the effective tidal radius for 
Monte Carlo simulations is rt^^f = rt/xud and it is smaller than 
rt. Therefore, for Monte Carlo simulations, a system was slightly 
too concentrated (as measured by the ratio between the tidal radius 
and the half-mass radius) compared to N-body simulations; (ii) the 
escape process is instantaneous. A star which energy is greater than 
the critical energy is promptly removed from the system; (iii) the 
coefficient xtid is an explicit func tion of A^. Its dependence wa s 
calibrated only for low N system ( IGiersz. Heggie. & Hurley|2008l) . 
so one can have some doubts about the rate of system evolution for 
large N. 

To overcome the drawbacks which are described above, we 
decide d to apply the theory described in iFokushige & Heggi3 
( I2OOOI) . According to this theory the time-scale for escape is given 
by 

7v{E~E,rur ' 
where uj is the angular velocity of a cluster around a parent galaxy 
and E is the energy of a star, ytid is a coefficient, which slightly 
depends on the system structure, and can be approximated by 0.38. 
The probability of escape in a time-step At of a star with energy 
greater than Ecru is given by 



Pa{At) = 1 - exp{-At/Q, (2) 

According to Figure 9 in lpokushige & Heggid ( 120001) the equation 
([2J matches the simulation results very poorly not only for the es- 
cape time scale but also for the overall shape of escape probability 
distribution. Indeed strictly it is known only to give an upper limit 
to the rate of escape, and they found empirically that the true rate 
of escape is smaller by about a factor of 10. So we can treat ytid in 
equation ([T) as a free parameter and adjust it by comparison with 
N-body simulations. 

To better represent the empirical shape of probability distri- 
bution, which can not be properly represent, even with appropriate 
choice of ytid, we decided to use approximation to equation ^ 
suggested by[Fokushige & Heggie (2000) for the probability distri- 
bution of escape times. For the purposes of the MOCCA code, this 
approach is strictly formal and does not have any physical meaning. 
The probability of escape in a time-step At is given by: 

P/(A^) = l-(l + a^)-^ (3) 

where a and b are coefficients which slightly depend on the 
structure of the system and i — ujAtE'^, where E = {E ~ 
Ecrit) /\Ecrit\- In the MOCCA code they are a free param- 
eters which are fitted against results of N-body simulations. 
As a first guess, values e qual to 3.0 and 0.8 can be adopted 
dFokushige & Heggiell2000l) . 

So, to mo del in MOCCA escape process according to 
IFokushige & Heg gie ( 2000) we have to adjust against N-body sim- 
ulations the free parameters: ytid or a and b. The big advantage of 
this approach is that the probability of escape does not explicitly 
depend on the number of stars in the system and that the escape 
process introduced into the MOCCA code generally follows the 
one in N-body systems. 

2.2 Probability of interactions 

In the hyperbolic 2-body approximation, the total cross section for 
interaction between a binary and a star, or another binary, with a 
pericentre distance less than rp,nax, is given by, 

2 2 ( , 2Gmi23 \ 
o- = 7rp = -nrp^ax 1 H ^ I , (4) 
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where p is the impact parameter, V is the relative velocity be- 
tween binary and star or other binary, mi23 is the mass of inter- 
acting objects. The second term in the brackets in equation ([4j 
describes the so coled gravitational focusing term. The larger the 
Tpmax, the larger the probability of interaction. In general, the en- 
ergy outcome from the interactions should not depend o n rj,max 
for r-p rnax ^ A, whcre A is a binary semi-major axis jHeggid 
For large enough rpmax there is a balance between positive 
and negative binary binding energy changes. So the tail of rpmax 
distribution is not important from the point of view of binary en- 
ergy generation. Of course r-pmax cannot be too large because in 
such a case part of interactions will be very similar to an ordinary 
relaxation process which is already covered by the Monte Carlo 
engine. On the other extreme, for small enough rpmax, practically 
all interactions will be resonant, hard and connected with strong 
energy generation. In this case substantial number of interactions 
with modest energy changes are missing. Of course, a larger rpmax 
means a larger impact parameter and larger number of interactions 
with the same time interval. Most of these interactions are very soft 
from the point of view of energy generation, but may lead to rel- 
atively large changes of binary eccentricity. This means that there 
is also larger probability for induced (during interaction process or 
shortly after it) mass exchange between binary components. This in 
turn can lead together with stellar/binary evolution to a larger rate 
of formation of "peculiar" objects (e.g. BSS) or a larger rate of bi- 
nary merger or disruption. Indeed this is the case as can be seen in 
Fig|7] Therefore, the determination of rpmax is not important from 
the point of view of total binary energy generation and evolution 
of cluster global parameters (e.g. total mass or core and half-mass 
radii), but is important from the point of view of the formation of 
many kinds of "peculiar" objects which are formed in direct inter- 
action between stellar dynamics and stellar evolution. An interest- 
ing influence of rpmax on the spatial distributions of some binary 
parameters is discussed in Sec l3.3.3l below. 

To roughly assess the range of values of rpmax we will com- 
pare the cross section given in equation (|4j with integrated differ- 
ential cross sections over all possible binding energy changes for 
resonant interactions according to Heggie formulae (Spitzer 198^, 
eq. 6-23), for flybys and resonant interactions according to Spitzer 



_ bys 

formulae i SDitzeJl987l. eq. 6-27) and for hard binary-binary inter- 
actions ( lGaoetalJl99ll. eq. 2.7). 

2AHmim2 



A 



= < 



7v37r 771127713 

hAs'Kmi_m2 



16\/37r7/7l2?7l3 

104 
42^ 



for Heggie 
for Spitzer (5) 
for Gao — equal mass case. 



where mi and 7n2 are masses of binary components, mi2 = 
mi + m2, ms is the mass of an approaching star, A is binary semi - 
major axis. As and Ah are coefficients defined in ISpitzeJ l fl987h . 
For an equal mass case the ratio rpmax/ A is equal to: 1.1, 3.4 and 
0.8 for the Heggie, Spitzer and Gao cases, respectively. The aver- 
age values of rpmax/ A obtained from the real simulations (com- 
puted for every interaction) according to equation (O are close to 
the above values. 



3 N-BODY - MOCCA COMPARISON RESULTS 

The detailed comparison between results of N-body and MOCCA 
simulations will proceed in three steps. First, the best values of 
rpmax, ytid, a and b will be chosen by the comparison of the total 



mass, the half mass radius, the core radius, the binary number and 
BSS number. Second, the results presented in^Baumgardt (20()l|) 
for small A'^ and single mass systems will be checked - the scaling 
of the half time with A*' and the evolution of the number of poten- 
tial escapers. Third, the detailed comparison for binary spatial and 
energy distributions, binary binding energy, evolution of some of 
Lagrangian radii, average masses in different parts of system and 
etc. will be done. The detailed comparison between N-body and 
MOCCA simulation results will help the reader to independently 
asses from our conclusions how well MOCCA can follow N-body 
and how reliable a code it is. 

The data for the comparison with MOCCA simul ations comes 
from simulations done by Jarrod Hurle y for M67 jHurlev et al.l 



l2005h with N = 24000, for NGC6397 l lHurlev et al]]2008l) with 
N = 100000 and for A^ = 200000 (private communication). 
The initial conditions for the N-body simulations are summarised 
in Tab^ 

To minimise the statistical fluctuations connected with the 
generation of the initial models the initial positions, velocities, 
masses and binary eccentricities and semi major axes for the 
MOCCA simulations are taken directly from the N-body simula- 
tions. They are provided as an input files. The observed in MOCCA 
simulations statistical fluctuations are connected only with different 
sequences of interactions and movements of objects in the systems. 
Positions, velocities and objects chosen for interactions are picked 
up according to the Monte Carlo technique. 

3.1 Determination of the Free Parameters 

The probability of escape given in equations ([2]and[3]( depends, 
among other factors, on the time. At. According to the Monte Carlo 
approach, t o scale ti me step from the Monte Carlo units to the N- 
body units (Heggie & Mathiieuiil986) one needs to use the follow- 
ing equation. 



AtNB = At 



MC 



N{t) 



(6) 



where N{t) is an actual number of objects in the system and 7 is 
the coefficien t in the Coulomb logarithm, for multi mass system 
equal to 0.02 l lGiersz. Heggie. & Hurlevll2008h . 

As can be seen from Fig.Q the scaling from Monte Carlo to 
N-body time units (according to equation {6}) gives inconsistent re- 
sults with N-body simulations - evolution of the total mass is too 
slow. This means that the escape rate in MOCCA simulations is 
considerably slower than in N-body. The same is true for other N 
N-body simulations. It has to be stressed that the both models can 
not be brought into agreement by appropriate choosing of any free 
parameters. Surprisingly, if we compute the escape probability by 
scaling the time step by a formula like equation ([6]l, but using the 
initial number of stars A'^(O) instead of the current value N{t), then 
it is found that a good match is obtained with N-body results, for all 
values of N that we have checked. The reason for this behaviour is 
u nclear There are at least two possible explanations: (1) according 
to lFokushige & Heggiel ( I2OO0I) the free parameters in equations (|2] 
and[3) depend on the structure of the system. For more concentrated 
King models the coefficients are larger During cluster evolution the 
structure of the system is changing (core - halo structure is devel- 
oping), so one can expect that the free coefficients can depend on 
time, (ii) The shape of the cluster around of the tidal radius depends 
on the energy of a star. This energy in the reference frame fixed to 
the cluster centre moving around a host galaxy with the angular 
velocity uj is given by: 
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N = 100000 Rp=2a-a 



N = 24000 Rp=2a-a 




Figure 1. The evolution of the system total mass as a function of time 
for different MOCCA models and N-body simulations. MOCCA-nt means 
scaling between Monte Carlo and N-body time units according to equation 
(6) using the current value of N{t), MOCCA-ntO means that, for the pur- 
poses of determining the escape probability only, scaling between Monte 
Carlo and N-body time units uses the initial number of objects in the sys- 
tem, N{0), in place of N{t) in|6] MO CCA-xtid=1.27 means results fo r 
scaling according to prescription given in lGiersz, Heggie, & HurlevlfcOOSi) , 
and N-body means N-body results. 



E : 



I ^ I 1 2. 2 



3x 



for N-body 
for MOCCA 



(7) 



where (j) is the potential v is the velocity dispersion, and x and z 
are coordinates with origin at the cluster centre. The last term for 
N-body depends on the centrifugal and tidal forces. In the case of 
the MOCCA code this term is not present, so the star energy is not 
exactly comparable in the both methods. The MOCCA system has 
spherically symmetrical shape instead of "oblate". This could lead 
to some differences in the time to escape between MOCCA and N- 
body code, which we attempt by overcame by choice of yud or a 
and b. 

The escape criterion used in the o ld version of the Monte Carlo 
code l lGiersz. Heggie, & HurlevlboOSh gives reasonable agreement 
with N-body results (despite it was calibrated only for low A^), 
although the rate of evolution is systematically to fast. This as- 
sures that the models of real star clusters computed with the 
old Monte Carlo code are relevant. Nevertheless, the MOCCA 
code with the new descrip tion of the escape process (based on 
iFokushige & Heggi3 dlOOOl) ') gives more consistent results with N- 
body, not only with respect to the evolution of the global parame- 
ters, but also with respect to the detailed properties of binary dis- 
tributions. What is also important, it is A'^ independent and can be 
safely used for any A*'. 

In the rest of the paper we use A(0) in the escape procedure 
instead of N{t) for scaling time from Monte Carlo to N-body units. 
For all other processes the scaling given in equation ^ was used. 

To determine the free parameters described above (ynd or a 
and b) we have to run several simulations with a different number 
of stars and different values for the parameters, and then compared 
results with N-body simulations. The best values were chosen "by 
eye", we did not attempt to asses how accurate and how unique they 
are (but see comments about statistical fluctuations below). 
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Figure 2. The evolution of the system total mass as a function of time for 
MOCCA models with different coefficients a and b and N-body simulations 
for N = 24000. 



N = 24000 Rp=2a-E 




2000 3000 

Time (Myr) 



Figure 3. The evolution of the system total mass as a function of time for 
MOCCA models with different coefficient ya^ and N-body simulations for 
TV = 24000. 



The results are presented in Figs|2l[3]and|4] As one can see, 
the dependence on yud is stronger than on a and b. The "best" 
values are yud ~ 4.0 and a — 3.0, 6 = 0.7. The values 
for a and b are close to the values given in IFokushige & Heggid 
in Table 1 there, ytid i s about ten times larger than given 



Fokushige & Heggid ( I2OOOI) (eq. (9) there), but taking into ac- 



count that the time scale given in equation ^ is too short (by about 
1 dex) the value 4.0 is ver y close to the empirical value which 
IFokushige & Heggid 1 I2OOOI) found. Only for the simulations with 
A^ = 24000 the evolution rate is slightly too slow for MOCCA, 
for time larger than about 2 Gyr, and it seems that yud = 3.0 
is a better choice than ynd = 4.0. For other A'' (100000 and 
200000) MOCCA simulations follow the N-body results very well 
and choice of yud ~ 4.0 gives better agreement with N-body. 

To asses the influence of statistical fluctuations on the obtained 
results, three simulations with exactly the same initial conditions 
but with different sequences of random numbers were run. The re- 
sults are given in Fig. |5]for N = 24000. The fluctuations for this 
model are largest, but clearly smaller than the difference connected 
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Figure 4. The evolution of the system total mass as a function of time for 
MOCCA models with different coefficients a and b and N-body simulations 
for N = 200000. 
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Figure 5. The evolution of the system total mass as a function of time 
for MOCCA models with a = 3.0 and b = 0.7 and different ran- 
dom sequences (iseed=10, iseed=20, iseed=30) and N-body simulations for 
N = 24000. 



with different values of the free parameters. For larger A'^ the fluc- 
tuations are practically negligible. For other global quantities the 
fluctuations are similar like for the total mass. They are small and 
presented in the paper results for a single simulations are represen- 
tative. 

Having determined a, b or yud we now turn to find out the 
best value for r-pmax- As it was argued above the value of r-pmax 
will have big impact on the number and distribution of binaries and 
BSS in the system. In the MOCCA code BSS is defined exactly 
as it was done in Hurley's N-body simulations. A main sequence 
star is identified as BSS when its mass is greater than 1.027l/to, 
where A'lto is the turn off mass. As can be seen from Figs[6]and|7] 
the requirements set by the N-body simulation for the number of 
BSS and binaries are rather contradictory from the point of view of 
MOCCA results. 

To get the best agreement for the evolution of the total number 
of binaries, MOCCA needs as large as possible r-pmax- rpmax = 
27 A seems to be a good choice. On the other hand to match N-body 
results for the number of BSS MOCCA prefers modest values of 
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Figure 6. The evolution of the number of binaries in the system as a func- 
tion of time for MOCCA models with different rpmax and N-body simu- 
lations for N = 100000. e.g. Rp=27a-9a means VpYnax 

= 27 A for three- 
body interactions and rpmax = 9^ for four-body interactions. 
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Figure 7. The evolution of the number of BSS in the system as a function 



of time for MOCCA models with different rp„ 
for TV = 100000. e.g. Rp=27a-9a means rp„ 
interactions and r„m 



■ and N-body simulations 
; = 27^ for three-body 



9A for four-body interactions. 



Tpmax - equal to about A. That conclusion is also true for = 24000 
and N = 200000. It seems that a reasonable compromise between 
the evolution of the total number of binaries and BBS is given for 
the rpmax suggested by the theoretical considerations given above 
in Sec |2.2| namely for three-body interactions rpmax ~ 2A and 
for four-body interactions rpmax = A. One can thing that this con- 
tradictory requirements set by the number of BSS and binaries can 
be explained by the different definitions of binary in N-body and 
MOCCA simulations. In the N-body results presented here, binary 
is identified with a regularised binary (rather hard), so called KS 
binary. In MOCCA we follow all binaries, even very soft ones. So, 
generally in MOCCA simulations we should expect larger number 
of binaries than in N-body. To quantify this we checked the number 
of non KS binaries in N-body simulations. This number was rather 
small (at most about 30) and could explain observed differences 
only partially. 
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It is worth to note that, as one can expect, for the MOCCA 
code the evolution of the Lagrangian radii and the total mass do not 
depend on the value of rpmax - the total energy generation in three- 
and four-body interactions does not depend on rpmax provided 
that rpmax is not too small or too large. The situation is differ- 
ent when the FewBody integrations are switch off in MOCCA and 
only cross sections are used for energy generation in binary inter- 
actions (MOCCA-NoFB). Then there is a very strong dependence 
on rpmax for the Lagrangian radii evolution. Larger rpmax means 
larger probability for interactions. Each interaction generates on av- 
erage the same amount of energy (according to the adopted cross 
section), so for larger rpmax more energy is generated by binaries, 
and the system expands faster than in the case of MOCCA sim- 
ulations (with FewBody integrator). It seems that the best values 
of rpmax for MOCCA-NoFB are 0.5A and 0.25A, for three- and 
four-body interactions, respectively. 

In the remainder of the paper the best values of the free pa- 
rameters given in this section yud = 4.0, a — 3.0, b = 0.7, 
rpmax = 2A and A for three- and four-body interactions, respec- 
tively will be used for comparison of different system and binary 
properties obtained in MOCCA and N-body simulations. 



(t,h) 



MG 

0.61 ''rl>'^' 



Figure 8. The dependence of the half mass time on the initial total number 
of stai's A'^. Green line the standard linear scaling with the half mass re- 
laxation time, t^h, blue line scaling predicted bv lBaumgardJ j200ll) , i^^*, 
purple line scaling obtained in this paper, t^^^ , pal blue line scaling with A'^ 
in this paper, N^^'^"^. Points with error bars (Scr) are the results of MOCCA 
simulations. 



3.2 Half mass time and potential escapers 

Heaving chosen the free parameters a, b and yud, which determine 
the rate of escape from the tidally limited cluster we can now check 
the results of si mulations for a tida lly limited and single mass sys- 
tem presented in lBaumgardtlj200lb . Namely: the A'^ dependence of 
the half-time, time when system contains half of its initial mass, 
and the evolution of the number of potenti al escapers for di ffer- 
ent A'^. According to the model presented bv lBaumgardll ( 1200 ih the 
half-time should scale as the relaxation time to the power 3/4 in- 
stead of the linear scaling with the relaxation time predicted by the 
standard theory. On Fig[8]is show the half-mass time as a function 
of N for MOCCA simulations from N = AK lo N = 256K. It 
also shows different scaling laws fitted to the simulation data. As 
we can see for A'^ up to 16 A" (the range for models presented in 
iBaumgardj l l200lh ) the scaling proportional to t^^* is acceptable 
within error bars (Sa - cr is estimated from 10, 7, 5 and 3 simula- 
tions with the same initial conditions but with different random se- 
quences for N = 4fc, N = 16k, N — 64k and N = 256fc, respec- 
tively). For larger N the scaling is less steep and the whole range 
of A*' can be reasonably well fitted with scaling 4°^^ or A"'^*. 
The reason for the discrepancies between results predicted by the 
theory scaling and those obtained from the simulations may be at- 
tributed to the facts that the star energy in MOCCA does not contain 
terms attributed to the non inertial reference frame a nd tidal force. 
The sc alin g A^'^-^'^ is close to the resul ts obtained bv lLamers et al.l 
( l2005l) and lGieles & Baumgardll | [200§) . Ar° ^ 

Fig|9]shows the evolution of the potential escaper fraction with 
time for different initial number of stars. The initial model setup 
is responsible for the 15 per cent potential escapers at the begin- 
ning. The cluster starts with primordial escapers, because the es- 
cape energy, Ecru, is lower than the edge potential of the initial 
King model. The different behaviour of the number of the poten- 
tial escaper for different A at the beginnin g is in agreement with 
the results presented bv lBaumgardJ hoOll) (his Fig. 11). The ini- 
tial increase of the potential escapers probably indicates the phase 
when cluster evolve towards equilibrium after removal in a short 
time substantial amount of mass feaumgardt 200 li) . This increase 
is largest for small N systems, for largest N is barely visible. The 





0.6 0.4 0.2 
N/N(0) 

Figure 9. Potential escaper fraction as a function of stars bound to the clus- 
ter for N = AK, N = IQK, N = 64K and N = 256X, clockwise. 



number of potential escapers decrease with time until the core col- 
lapse, when it starts to rise again. The c omparison of Fig |9] for 
A' = 4fc and A' = 16k with Figure 11 in lBaumgard3 j200lh sug- 
gests that the core collapse is delayed for MO CCA. The detailed 
inspection with Figure 8 in iBaumgardJ ( 1200 ih shows that the de- 
lay is rather modest not larger than 5—10 per cent. The reasons 
for this can be connected with the fact that the rpmax = A cho- 
sen for simulations was to small and the probability for binary for- 
mation in the three body interactions can be slightly too small in 
MOCCA. Indeed, larger rpmax = 27 A and larger formation prob- 
ability both brings MOCCA results to much better agreement with 
N-body, but still core collapse in MOCCA is slightly delayed. This 
suggests that there are other factors which could be responsible for 
this disagreement, e.g. the 7 coefficient in the Coulomb logarithm, 
or the rate (in N-body and MOCCA) in which small A clusters 



© 2002 RAS, MNRAS OOO.fTlflSl 



8 M. Giersz, D.C. Heggie, J.R. Hurley & A. Hypki 



can regain equilibrium after the initial substantial mass loss. The 
observed disagreement does not influence the results for the half 
mass time and the evolution of the potential number of escapers 
before the collapse time. We can conclude that the MOCCA code 
reproduce reasonably well results presented bv lBaumgardj bOOlh 
for small N-body simulations. 



3.3 Results of comparison 

The comparison between MOCCA and N-body results will proceed 
in three steps. First, the evolution of the global parameters (total 
mass, Lagrangian radii, core radius) will be compared. Second, the 
evolution of the binary properties (number, energy, mass and num- 
ber distributions) will be checked. Third, properties of the pecu- 
liar objects like BSS and black holes (BH) will be compared. Most 
of the figures presented below, will additionally display the result 
of MOCCA simulations with a switched off FewBody integrator, 
and a switched on interaction cross sections (MOCCA-NoFB). This 
will help the reader asses how well the simplified MOCCA-NoFB 
(very similar to the old version of the Monte Carlo code, which was 
previously used to successfully simulate evolution of real star clus- 
ters) can follow N-body results and for which cluster properties it 
is enough to use the much faster and simplified code. 



3.3.1 Global parameters 

The comparison between N-body and MOCCA results was par- 
tially discussed already in Sec B.ll for the total mass evolution. It 
was shown that the agreement between both techniques is very 
good, only for A'^ — 24000 it was barely acceptab le (see Fig|3). The 
evolution of the core radius (defined according to lCasertano & Hti3 
( Il985h ) for TV = 200000 is shown in the Fig[TO] The agreement 
between MOCCA and N-body is very good. As one can expect, 
MOCCA-NoFB gives a slightly too large core radius. This, as it 
was explained in Sec |3.1| is connected with the overestimation of 
the binary energy generation in the cross section regime for a large 
enough rpmax. The large fluctuations in the core radius visible in 
the figure for N-body and MOCCA simulations are connected with 
the movement of a massive BH or BH-BH binary in the system. 
The mass of the BH objects is about 30 — 50Mq. Its movement 
in the system is connected with kicks acquired in interactions. If 
the massive object is present in the core the core radius is smaller, 
than when it is in the halo. The mass of the massive object is com- 
parable to the core mass. When all the most massive objects are 
removed from the system, because of strong interactions with other 
massive binaries or stars, the evolution of the core radius is again 
"smooth". The evolution of the core radius is generally similar for 
the other models (but without a large fluctuations) with exception 
for = 100000 when N-body results are systematically slightly 
below the MOCCA (see ¥ig.^. 

The evolution of the half mass radius for A*' — 100000 model 
is shown on the Fig[TT] The evolution of all models from the very 
beginning is very similar, although the N-body is slightly below 
both MOCCA models. The differences starts to build up around 
the core collapse time (bump), which seems to be around 17 Gyr for 
N-body. For MOCCA models the core collapse time is around 20 
Gyr. The bump in the half mass radius is also visible for MOCCA, 
but less pronounced. For N = 24000 and N = 200000 the half 
mass radius for MOCCA models is systematically slightly below 
the N-body models. The differences start to build up from the be- 
ginning, and are biggest around the time when the stellar evolution 
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Figure 10. Evolution of the core radius (detined according to 
ICasertano & Hull h985l) ) for N = 200000. Red line - MOCCA, green 
line - MOCCA-NoFB and blue line - N-body. 
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Figure 11. Evolution of the half mass radius for N = 100000. Red line - 
MOCCA, green line - MOCCA-NoFB and blue line - N-body. 



stops to be the dominant process of cluster expansion (time when 
the indirect heating connected with the stellar mass loss becomes 
smaller than the heating connected with the binary energy genera- 
tion). Then the evolution of the half mass radius starts to converge 
for both N-body and MOCCA models. The same behaviour can be 
observed for other Lagrangian radii (1% and 10%). The reason for 
such a behaviour is unclear and can be attributed to different mu- 
tual strengths of the physical processes which operate during the 
different phases of cluster evolution. Faster mass segregation in N- 
body than in MOCCA can generate more extended cluster. Stellar 
evolution responsible for the loss of stellar mass can substantially 
blow up the cluster, particularly at the initial phases of evolution. 
Both N-body and M OCCA models rela y on the same stellar evolu- 
tion p rescription (H urley. Pols. & Tout 2000; Hurley. Tout. & Polsl 
I2OO2I) . so we cannot expect that the amount of mass loss is differ- 
ent in both models. However, the mass segregation acting together 
with the stellar mass loss can substantially amplify the expansion 
effect. If the most massive stars lose their envelopes when they are 
already mass segregated the effect on cluster expansion is largest. 
Finally, the larger binary energy generation in N-body simulations 
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Figure 12. Evolution of the average mass inside 50% Lagrangian radius for 
N = 200000. Red line - MOCCA, green line - MOCCA-NoFB and blue 
line - N-body. 
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Figure 13. Evolution of the total binding energy of binaries for = 
200000. Red line - MOCCA, green line - MOCCA-NoFB and blue line 
- N-body. 



(see Fig[T3]- increase of the total binary binding energy) can be re- 
sponsible for faster Lagrangian radii expansion. The evolution of 
the average mass inside 50% Lagrangian radius and evolution of the 
total binary binding energy for A'^ = 200000 are shown on Figsll2l 
and [T3] These figures are representative for all models and La- 
grangian radii. Indeed, it seems that the mass segregation is slightly 
stronger in N-body than in MOCCA. This suggests that the larger 
mass segregation in N-body model can be responsible, at least par- 
tially, for the slightly discrepant evolution of the Lagrangian radii. 
The evolution of the total binary binding energy, is from the very 
beginning rather similar, until later time (about 4 Gyr), when more 
bound binaries are formed in N-body. The final formation of very 
hard binary is visible in MOCCA and N-body (see discussion in 
Sec irm , but not for MOCCA-NoFB. When binary is removed 
from the system, because of interactions, all models again are sim- 
ilar. It is worth to note that the results of simulations for MOCCA 
and MOCCA-NoFB are very similar until the late phases of evolu- 
tion. 



3. 3. 2 Binary properties 

The evolution of the number of binaries was already presented in 
Sec |3.1 [ during discussion about the determination of the free model 
parameters. We know that the evolution of the total number of bi- 
naries is slightly too slow for MOCCA. The difference starts to 
buildup at the time when stellar evolution becomes less and less 
important (see Fig[6l(. The average binary mass and the binding en- 
ergy distributions will be discussed for model A'^ — 200000. The 
results for this model are representative for other models and ad- 
ditionally shows buildup of the average binary mass and binding 
energy for massive binaries. This buildup is possible only for the 
model N = 200000 for which the distribution of supernova (SN) 
kicks is uniform between and 100 km/s and lets larger, than in 
other models (Maxwellian distribution with a = 190km /s for SN 
kicks), fraction of BHs to be bound to the system. 

The evolution of the average binary mass for different regions 
of the system is shown in Figs[T4l[T5]and[T6] The agreement be- 
tween N-body and MOCCA results is very good in all regions. This 
is despite the fact, that in N-body there is a smaller number of bi- 
naries than in MOCCA and Lagrangian radii are slightly different. 
It seems that the average binary mass and its distribution does not 
depend on the number of binaries, which is probably a result of ex- 
actly the same binary initial conditions for the both models and sim- 
ilar mass spectrum of removed or destroyed binaries in the models. 
For the region inside the core, the buildup of the mass of binaries is 
clearly visible. The big fluctuations are connected with the move- 
ment of the massive binaries, which because of hard interactions 
with other stars are kicked out from the core and then because of 
mass segregation sink again to the centre. Finally, they are kicked 
out of the system and the average mass of binaries starts to become 
less chaotic and its changes are rather small. The drop and then in- 
crease of the average mass around I4.5Gyr is connected with the 
core collapse. The high central density makes binary-binary inter- 
actions very effective and substantial number of relatively wide and 
massive binaries are destroyed causing drop of the binary average 
mass. Looking on the average masses in the different cluster re- 
gions the mass segregation is clearly visible. In the centre the mass 
is about twice as large as in the halo. For the A'^ — 200000 model 
the difference between MOCCA and MOCCA-NoFB is very small, 
except in the core when the increase of the average mass is less 
pronounced for MOCCA-NoFB. The sharp increase of the aver- 
age mass, close to the cluster dissolution time, is connected with 
the fact that only the most massive binaries are able to stay in the 
system, less massive ones are successively removed. 

For other N models the behaviour of the average binary mass 
is similar to the one described above, except that in the core there is 
no large fluctuation (there are no very massive binaries there) and 
for late phases of cluster evolution small discrepancies between N- 
body and MOCCA start to show up. 

The evolution of the average binary binding energy for dif- 
ferent regions of the system is shown in Figs[T7][T8]and[T9] The 
best agreement between N-body and MOCCA results is for the 
core region. The buildup of the binary binding energy connected 
with the formation of massive BH-BH binary is clearly visible. It 
has to be stressed that the increase of the binary binding energy is 
not connected with the core collapse and core bounce. It is purely 
connected with the formation of a very massive BH-BH binary and 
then in the increase of its binding energy in interactions with stars 
and other binaries. The binary is finally removed from the system 
and then the average binary binding energy suddenly drops. For 
other cluster regions the agreement between N-body and MOCCA 
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Figure 14. Evolution of the average binary mass between 50% Lagrangian 
radius and the tidal radius for = 200000. Red line - MOCCA, green line 
- MOCCA-NoFB and blue line - N-body. 
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Figure 16. Evolution of the average binary mass between 50% Lagrangian 
radius and the core radius for A'^ = 200000. Red line - MOCCA, green line 
- MOCCA-NoFB and blue line - N-body. 



I^OCCA - Rp=2a-a 
I^OCCA-NoFB - Rp=2a-a 
N-body 




8000 10000 
Time (Myr) 



Figure 15. Evolution of the average binary mass between 10% and 50% La- 
grangian radii for TV = 200000. Red line - MOCCA, green line - MOCCA- 
NoFB and blue line - N-body. 
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Figure 17. Evolution of the average binary binding energy between 50% 
Lagrangian radius and the tidal radius for N = 200000. Red line - 
MOCCA, green line - MOCCA-NoFB and blue Hne - N-body. 



is less satisfactory. The average binary binding energy for N-body 
is systematically larger than for MOCCA. The differences start to 
buildup just after the time when the stellar evolution stops to dom- 
inate the cluster global evolution. It seems that in N-body simula- 
tions harder binaries can be produced in the core before they are 
kicked out to the outer parts of the system. Maybe this is connected 
with the fact that in MOCCA code only binaries are allowed, the 
higher hierarchies are not allowed. Triples and quadrupoles are ar- 
tificially disrupted to binaries and single stars. It is well known that 
in N-body simula tions substantia l number of triples and quadruples 
are formed (e . g. lMikkolalll984 iMcMillan. Hut & Makinclll99TI ; 



iHeggie & Hu3 llOOj iHurlev et alj l20oT and reference therein) 
They can interact with other object in the system and produce, on 
average slightly harder binaries. Differently than for the average 
binary mass the average binary binding energy for MOCCA-NoFB 
is clearly smaller than for MOCCA and N-body. In MOCCA-NoFB 
simulations Heggie's cross section was used (see Sec |2.2| ). There- 
fore, on average the binary binding energy is changing by 40% in 
every interaction. It seems that such large energy changes are too 



big to produce a substantial population of binaries with energies 
smaller than the energy needed for binary escape from the system 
and larger than the minimum binding energy needed to generate bi- 
nary escape in next interactions. Smaller energy changes will allow 
binaries to populate such a binding energy region. For the MOCCA 
simulations with rpmax ~ '2A for binary-single interactions, the 
average change of binding energy is about 18%. It seems, from the 
point of view of under abundance of binaries with large enough 
energies in the outer parts of the system, that this number is still 
too large. Indeed, larger rpmax gives much better agreement with 
the N-body results (the larg the smaller average change 

of binary binding energy), but is still not perfect. We should keep 
in mind that forcing too large r-pmax to get smaller binding en- 
ergy changes in interactions we create a much larger number of 
BSS comparable to N-body simulations (see Sec l2.2l and Fig.|7]l. In 
this paper we decided not to use extremely large values of rpmax, 
putting more emphasis on the number of BSS than on the distri- 
butions of the average binary binding energy. The dependence of 
the number of BSS on rpmax is much stronger than for the bi- 
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Figure 18. Evolution of the average binary binding energy between 10% 
and 50% Lagrangian radii for = 200000. Red line - MOCCA, green 
line - MOCCA-NoFB and blue line - N-body. 
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Figure 20. The evolution of the number of blue stragglers in the system as 
a function of time for MOCCA models with different Tpmax and N-body 
simulations for N = 24000. 
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Figure 19. Evolution of the average binary binding energy between 50% 
Lagrangian radius and the core radius for = 200000. Red line - 
MOCCA, green line - MOCCA-NoFB and blue line - N-body. 



nary binding energy. But if one is more interested in binary bind- 
ing energy distribution, much larger Vpmax, say 9 A or even higher, 
should be used. It is worth to stress that MOCCA-NoFB cannot re- 
produce the binary binding energy distribution observed in N-body 
and MOCCA. 

For other N-body models the behaviour of the average binary 
binding energy is similar to the one described above, except that in 
the core there is no large buildup of very hard binaries. For A'^ = 
24000 the discrepancies between N-body and MOCCA are largest. 
This is probably connected with the fact that this model contains as 
much as 50% binaries (10, 20 times more than other models), and 
small differences in binary energy generation may produce larger 
discrepancies in properties of binary distributions. 



3.3.3 Black holes and blue stragglers 

N-body simulations used for the comparison with MOCCA pro- 
vided additionally information about global system parameters and 
detailed data about binaries and about BHs and BSS. Some data on 



BSS was already used in Sec B.ll to put some constrains on rpmax. 
In this section we will try to asses how well MOCCA code can re- 
produce the N-body results in respect to behaviour of "peculiar" ob- 
jects, which are rare and are formed in interaction channels, which 
are not very probable. 

In Fig|20] evolution of BSS for different rpmax for A'^ — 
24000 is shown. The dependence of the number of BSS on Vpmax is 
different than mN = 100000 and iV = 200000. There is no sharp 



increase of the maximum number of BSS for r„ 



> 2,A-A (see 



FiglT}. For N = 24000 there is a very weak, if any, dependence of 
the number of BSS on rpmax- The fluctuations of the number of 
BSS connected with the stochasticity of the system and MOCCA 
model are comparable in the size with differences related to rpmax 

- about 5 BSS. The differences observed on Figs|7]and|20]are rather 
surprising and point in the direction of strong dependence on the bi- 
nary fraction. It was 50% m N = 24000, at least ten times more 
than in the other models. This result contradicts the naive think- 
ing that a large number of binaries should amplify the effect of 
large rpmax- Maybe, a large binary fraction makes binary-binary 
interactions more probable and more binaries are destroyed or bi- 
nary properties are changed in such a way that formation of BSS 
is less probable, somehow balancing the larger BSS formation for 
larger rpmax- Generally, MOCCA is able to reproduce the evolu- 
tion of the number of BSS given in N-body simulations provided a 
reasonable values of rpmax is chosen. MOCCA-NoFB gives much 
smaller number of BSS than other models. That is not surprising, 
because important for BSS formation channels are missing in the 
cross section approach. 

There is a long debate about BH kick velocities in SN ex- 
plosions. We tested three possible variants in the paper (guided by 
choices done in N-body si mulation s ): Firs t, Maxwellian distribu- 
tion with a = 190 km/s IPhinnevI l ll992l) . and Maxwellian dis- 
tribution with (7 = 190 km/s but with the maximum kick veloc- 
ity set to 50 km/s (reduction of the kick velocity to maximum 50 
km/s is the procedure used in Jarrod Hurley's N-body simulation) 

- kfallb — and kfallb = 2 in the figures, respectively. Sec- 
ond, uniform distribution with kick velocities between and 100 
km/s, and uniform distribution with kick velocities between and 
50 km/s - kfallb = and kfallb =2 in the figures, respectively 
(only for N = 200000 model). Third, Maxwellian distribution 
with a = 190 km/s but the kick velocity is finally modified by 
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N = 1 00000 Rp=2a-a 



MOCCA - kfallb=0 
MOCCA - kfallb=1 
MOCCA - kfallb=2 



N = 100000 Rp=2a-a 



MOCCA - kfallb=0 
MOCCA - kfallb=1 
MOCCA - kfallb=2 



Time (Myr) 

Figure 21. Evolution of the number of BHs in the system for = 100000. 
Red line - kfallb = 0, green line - kfallb = 1 and blue hne - kfallb = 2. 
For cases kfalb = and 2 practically all BHs escaped from the system just 
after their formation. 



Time (Myr) 

Figure 22. Evolution of the number of BH-BH binaries in the system for 
N = 100000. Red line - kfallb = 0, green line - kfallb = 1 and blue 
line - kfallb = 2. For cases kfalb = and 2 there is no BH-BH binaries 
because no BHs left after SN kicks. 



the amount of mass which falls back on BH during SN explosion 
dBelczvnski. Kalogera & Buhkll2002h - kfallb = 1 in the figures. 
The amount of mass fall back depends on the core mass of a star, 
just before the supernova explosion. For the core mass greater than 
7.6M0 the kick velocity is km/s. 

The discussion of the influence of the SN kick velocity defini- 
tion on the system evolution and on properties of BH and BH-BH 
binaries will be presented for the A'^ — 100000 model. For other 
models the different kick velocities do not have strong influence on 
the system parameters. For — 24000 the escape velocity is very 
small, so even for kfallb = 1 case only a few BH will stay in the 
system after SN kicks and they will be quickly removed by inter- 
actions from the system. For A'' = 200000 the SN kick velocities 
for kfallb = and kfallb — 2 cases were already relatively small 
and a substantial fraction of BHs were left in the system. There- 
fore in the case kfallb — 1 the BHs fraction did not change sub- 
stantially and the system evolution was not influenced strongly. For 
A'^ = 100000 the situation is different. Only in the case kfallb = 1 
substantial number of BHs left in the system strongly influencing 
its evolution. 

On Figs|2T]and|22]the number of BHs and BH-BH binaries are 
shown for A'^ = 100000 and for different description of the SN kick 
velocities. The striking feature of these figures is the fact that there 
is no BHs and BH-BH binaries left in the system after SN kicks for 
case kick — and even for kick = 2. Kick velocities are too large 
and all BHs escape immediately from the system. Only in the case 
kfallb = 1 kick velocities are small enough to keep a substantial 
number of BHs in the system. The most massive BHs have masses 
larger than 25Af q , so they very quickly form in three-body inter- 
actions massive binaries with mass about SOAf© . The probability 
of binary formation in the three-b ody interactions is a very strong 
function of masses ( lHeggielll975l) . As can be seen from Fig |22l the 
number of BH binaries present at the same time in the system can 
be as large as 5. They are the most massive objects in the system, 
but this does not mean that all those binaries are in the same time 
in the core. They strongly interact with other stars/BHs or binaries 
and acquire kicks which move them outside the core. When they 
become hard enough they are one by one removed from the system 
by a final strong interaction. In the course of interactions each mas- 
sive BH-BH binary can remove several stars, therein a few BHs. 



N = 100000 Rp=2a-a 




5000 10000 15000 20000 



Time (lylyr) 

Figure 23. Evolution of the core radius for N = 100000. Red line - 
kfallb = 0, green line - kfallb = 1 and blue line - kfallb = 2. 

There are acting like a vacuum cleaner quickly reducing the num- 
ber of single BHs in the system, which is clearly visible in Figsl21l 
and|22] The probability of binary-single interaction depends on the 
mass of interacting objects. The larger the mass the larger the prob- 
ability and the larger number of hard interactions connected with 
substantial energy generation. The hardest and most massive bina- 
ries can produce much more energy in interactions than other bina- 
ries. Therefore, one can expect that such a large amount of energy 
is able to change the system structure showing up distinct features 
not visible in models with kfalb — and kfalb — 2. 

On Figs|23] and |24] are shown evolution of the core and the 
half mass radii for different kfallbs and N-body model for A'' — 
100000. We can see, that in the case when massive BHs and mas- 
sive BH-BH binaries are present in the system, the evolution of the 
core radius is very stochastic, characterized by large and fast fluctu- 
ations, on the time scale of an order of several dozen Myr. Also, the 
half mass radius in this case evolves faster and has larger values, 
which is presumably connected with stronger energy generation. 
The consequence of larger energy generation is much faster system 
evolution. The cluster lives shorter. Therefore, not only BH with 
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Figure 

kfallb 
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Time (Myr) 

24. Evolution of the half mass radius for A'^ = 100000. Red line ■ 
= 0, green line - kfallb = 1 and blue line - kfallb = 2. 



masses of several hundred AIq (IMBH) can influence the system 
structure, which can be detected by observations, but also presence 
in the system of several massive stellar mass BHs can leave obser- 
vational imprint on the cluster structure. That is a very interesting 
po ssibility worth to further check. Similar behaviour was reported 
m iHurlevI JioOTh for stochastic core rad ius evolution powered by 
single rather massive BH and in. Merritt e t al. ( 2004 ), Mackev e t al., 
( l2008h for strong core expansion of massive star cluster powered 
by large number of stellar mass BHs (in this case the core radius 
evolution is rather smooth because of a large number of BHs which 
forms bound subsystem in the central part of the system). 

The behaviour of the half mass and the core radii described 
above for N = 100000 is also visible for iV = 200000. The fact 
that in this case, for different kfallbs the number of BHs and BH- 
BH binaries are very similar, is interesting. This suggests that not 
the number of BH-BH binaries is crucial for the discussed above 
behaviour but the mass of BH-BH binary. For kfallb = 1 case 
more massive BHs can be formed than in other kfalbs because of 
the mass fallback during SN explosion. Therefore, the most mas- 
sive binary which can form in this case is substantially more mas- 
sive than in other cases. Such a binary can generate more energy 
in interactions and remove more stars and binaries from the sys- 
tem than a less massive binary. It is possible that in the case when 
two or more massive BH-BH binaries are formed in the system, 
they will interact between themselves and quickly harden. There is 
a non zero probability that before they escape, they will merge in 
collision interactions with "ordinary" stars/binaries and form more 
massive BH, which very quickly will form again binary. The se- 
quence can be repeated a few times forming a seed IMBH in the 
system. That is an interesting new road for the possible creation of 
IMBH in globular clusters. This possible scenario strongly relay on 
the SN kick velocity distribution for BHs, initial mass function for 
massive stars and initial cluster concentration. 

It is worth to stress that fast MOCCA-NoFB code can be used 
for projects which investigate the evolution of star clusters from the 
point of view their global properties. If one is interested in proper- 
ties of "peculiar" objects and their distributions should use slower 
MOCCA code. 

At the end of this section we would like to refer to the choice 
of the best value of rpmax- We argued that from the point of view of 
the number of BSS in the system the best choice of rpmax is about 



A, but from the point of view of the number of binaries or the dis- 
tribution of the binary binding energy larger values of about 9 A are 
more appropriate. The best compromise, taking into account statis- 
tical fluctuations, seems to be rpmax = 2A. This means that we put 
more emphasise on the limitations given by BSS than ones given by 
binaries. The reason for that is as follows. The BSS definition in N- 
body and MOCCA is exactly the same and the fluctuations in BSS 
number are smaller than the its value itself. So the sharp increase 
in the maximum number of BSS with larger rpmax (observed in 
FiglTJ is real and can not be attributed only to fluctuations. For the 
evolution of the number of binaries in the system and for the binary 
binding energy distributions the statistical fluctuations are practi- 
cally negligible, so we could expect that observed dependences on 
rpmax are real. But there are some doubts about the binary def- 
inition in N-body and MOCCA. In N-body simulation it is easy 
to count regularised binaries (KS binaries), but it is not quite so 
straightforward to find all non KS binaries. In MOCCA all bina- 
ries are directly followed. As we pointed out in Sec |3.1| however, 
the number of non KS binaries is rather small and can not be en- 
tirely responsible for the observed differences. For binary binding 
energy distribution evidences are rather indirect and the observed 
dependence on rpmax could be also, at least partially, attributed to 
unidentified yet physical processes or some systematical errors in 
MOCCA code. Therefore, we decided to opt for the compromise 
vale of rpmax than use larger values suggested by the total number 
of binaries and their distributions. We estimate that the error made 
in the total number of binaries by use of this compromise value is 
of order 5 per cent of the number of binaries 



4 CONCLUSIONS 

In this paper we have presented an advanced Monte Carlo code 
(MOCCA) for the evolution of rich star clusters, including most 
aspects of dynamical interactions involving binary and single stars, 
internal evolution of single and binary stars and complicated pro- 
cess of escape in the static tidal field. The direct integration of a 
few body pro blem was introduced on the base of FewBody code 
developed bv lFregeau et alj | |2004|) . The stell ar and binary internal 
evolu t ion was done according to BSE code ( iHurlev. Pols. & Toull 
I2OOOI : iHurlev. Tout. & Polsll2002h . The desc ription of the escape 
process was based on the theory presented bv lFokushige & Heggid 
I2OOOO. So, the MOCCA code is able to follow all channels of in- 
teraction up to binary-binary hierarchy including merging of stars, 
and the escape is not immediate any more, stars need time to find 
their way to escape from the Li and L2 Lagrangian points. The 
probability of escape and probability for interaction are character- 
ized by some free parameters which were adjusted by comparison 
of MOCCA and N-body simulation results for large A'^ systems, up 
toN = 200000. 

It was shown that the free parameters of the MOCCA code 
can be successfully calibrated against N-body simulations and that 
the free parameters do not depend on TV. MOCCA code is not only 
able to follow evolution of the cluster total mass, Lagrangian radii 
and the core radius, but also is able to reproduce in an reasonably 
accuracy distributions of binary parameters and number of BH-BH 
binaries and BSS. It also reproduces well the results obtained by 
iBaumgard^ 1 I2OOII) for single mass tidally limited systems for the 
half-mass time and evolution of potential escapers. The code is able 
to cope with very diverse systems, from single mass, isolated sys- 
tem without primordial binaries to multi mass, tidally limited sys- 
tem with large fraction of binaries. It was shown that simplified and 
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faster version of the MOCCA code (without direct FewBody inte- 
grator - MOCCA-NoFB) is a method of choice for projects which 
aim is to investigate the evolution of star clusters from the point of 
view their global properties. For other purposes, particularly when 
properties of "peculiar" objects and their distributions are in area 
of interests one should use slower MOCCA code. It is worth to 
note that MOCCA and MOCCA-NoFB simulations presented in 
this paper needs only about three and two hours, respectively, to be 
completed on a standard Opetron 2.4Ghz CPU. 

Despite these successes the MOCA code has still some known 
shortcomings, which we summarise here. 

(i) Higher-order multiples: It is widely argued that primordial 
triples and higher multiples should be incorporated into simulations 
along with primordial binaries. In any case, hi erarchical triple s 
form abundantly in binary-binary interactions ( iMikkolal Il984l ). 
Such higher-order multiples are ignored in the pres ent Monte-Carlo 
code, they are counted and artificially disrupted jHvpki & Giersj 

Hierarchical triples and higher-order multiples are planed to 
introduce in the new version of the code. 

(ii) Rotation: the Monte Carlo code is based on spherical sym- 
metry, and would require rather fundamental and very difficult re- 
construction in order to cope with cluster rotation. Rotation ac- 
celerates the rates of core co ll apse and mass segregation (e.g. 
iFiestas. Spurzem^ Kiml I2OO6I : I Kim et and references 
therein). In our models the absence of these rotational effects can 
be compensated by a modest alteration of the initial conditions. 

(iii) Static tide: the effect of tidal shocks have been extensively 
studied and it would be possible to add the effects as another pro- 
cess altering the energies and angular momenta of the stars in 
the simulations. The addition of tidal shocks will be more impor- 
tant when modelling Galactic globular clusters than open clusters, 
which usually are confined inside the Galactic disk. 

Despite these limitations, some of which are difficult to cure, 
the MOCCA code presented in this paper shows its potential power 
in simulations of real star clusters, from open clusters to rich glob- 
ular clusters. Monte Carlo models are feasible in a reasonable time 
(a few days) for globular clusters, which are, still for some time, too 
large for direct N-body simulations. The MOCA code is able to pro- 
vided data as detailed as N-body code can do. Only those two codes 
can provides such comprehensive information. Even when A'^-body 
simulations eventually become possible, Monte Carlo models will 
remain as a quicker way of exploring the parameter space for the 
large scale A^-body simulations. 
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